%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Jaynes-Cummings model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LM_Quantumfields_JcummingsDispersive
    global X0 Y0 Z0
    figure(1);          clf;
    
    o =2*pi*3;                                                              % coupling constant
    w =0;                                                                   % frequency
    d =2*pi*100;                                                            % detuning
    T =.2;                                                                  % time

    B0=[.5 .5]';                                                            % initial atomic state
    X=[0 ; 0 ; 0];                                                          % initial Bloch vector
    N =[1:30];                                                              % cut-off photon number
    n0=10;              a =sqrt(n0)*exp(0*1i*pi);                           % mean photon number
%     V0=sqrt(.5)*((n0==N)+((2*n0)==N));                                    % Fock state
    V0=exp(-n0/2)*a.^N./sqrt(cumprod(N));                                   % coherent state
%     V0=sqrt(.5)*(1+(-1).^N).*V0;                                            % cat state
%     V0=sqrt(n0.^N./ (n0+1).^(N+1));                                       % thermal state
    C0=B0*V0;           C0=C0/norm(C0);                                     % initial Jaynes-Cummings state

    x =pi*[0:.1:2]';
    X0=sin(x);          Y0=cos(x);          Z0=zeros(size(x));              % unit circles for Bloch sphere
    AA=n0*[-1:.02:1];   A =AA'*ones(size(AA))+1i*ones(size(AA'))*AA;        % complex plane for field states

    for j=1:200
        C = JCM_prop(C0,o,d,w,(j-1)*T,N);
        [X,Q,V]=JCM_calc(X,C,j,A,N);
%         W = JCM_wigner(C,j,A,N);
        W = A;
        JCM_plot(X,A,Q,W,V,T,N);
        pause(.01);
    end
end

function C = JCM_prop(C0,o,d,w,t,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hamiltonian and Schrödinger equation

    for n=N
        H =[(n-1/2)*w-d/2   o/2*sqrt(n)
            o/2*sqrt(n)     (n-1/2)*w+d/2	];

        C(:,n) = expm(-1i*H*t)*C0(:,n);
    end
end

function [X,Q,V] = JCM_calc(X,C,j,A,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculates Bloch vector, Q-function, and photon distribution


    X(:,j)=[sqrt(2)*real(sum(C(1,:).*C(2,:).''))                            % Bloch vector
            sqrt(2)*imag(sum(C(1,:).*C(2,:).''))
            sum(C(2,:).*C(2,:).''-C(1,:).*C(1,:).'') ];
%     norm(X(:,j))

    Q1=0;         Q2=0;                                                     % Q-function
    for n=N
        Q1=Q1+A.^n/sqrt(prod(1:n))*C(1,n);
        Q2=Q2+A.^n/sqrt(prod(1:n))*C(2,n);
    end
    Q = exp(-abs(A).^2).*(abs(Q1).^2+abs(Q2).^2);
  
    
    V = C(1,:).*C(1,:).''+C(2,:).*C(2,:).'';                                % photon distribution
end

function W = JCM_wigner(C,j,A,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% JCM in Wigner representation

    xA=real(A);         pA=imag(A);
    W =0;
    for m=N
        for n=N
            He_n_pA=2^(-n/2).*hermiteH(n,pA/sqrt(2));
            He_m_pA=2^(-m/2).*hermiteH(m,pA/sqrt(2));
            He_n_xA=2^(-n/2).*hermiteH(n,xA/sqrt(2));
            He_m_xA=2^(-m/2).*hermiteH(m,xA/sqrt(2));

            W = W+1/pi^2 * (C(1,m)'*C(1,n)+C(2,m)'*C(2,n)) ...
              .* (-1)./A * sqrt(pi/2/prod(1:n)/prod(1:m)) ...
              .* (-1)^(floor((m-n)/2)) ...
              .* (exp(-pA.^2/2).*He_n_pA.*He_m_pA+ ...
                  exp(-xA.^2/2).*He_n_xA.*He_m_xA);
        end
    end
end 

function Q = JCM_plot(X,A,Q,W,V,T,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot routine
    global X0 Y0 Z0

    subplot(241);plot3([X0 X0 Z0],[Y0 Z0 X0],[Z0 Y0 Y0],'b',X(1,:),X(2,:),X(3,:),'r')
    set(gca,'FontSize',16)
    xlabel('\fontsize{16}Re{\it\rho}_{12}')
    ylabel('\fontsize{16}Im{\it\rho}_{12}')
    zlabel('\fontsize{16}{\it\rho}_{22}-{\it\rho}_{11}')
    text(-2,-1,1,'\fontsize{16}(a)');
%     title('\fontsize{12}Bloch vector');
    axis tight
    view(30,50);
    subplot(242);plot([1:length(X(1,:))]*T,X(1,:),'b')
    set(gca,'FontSize',16)
    xlabel('\fontsize{16}{\itt}');	
    ylabel('\fontsize{16}Re{\it\rho}_{12}');	
    text(-10,.7,'\fontsize{16}(b)');
%     title('\fontsize{12}collaps and revival');
    axis tight;		
    subplot(243);bar(N,V,'r');		
    set(gca,'FontSize',16)
    xlabel('\fontsize{16}{\itn}');	
    ylabel('\fontsize{16}{p_n}');	
    text(-18,.15,'\fontsize{16}(c)');
%     title('\fontsize{12}photon distribution');
    grid;		
    subplot(244);mesh(real(A),imag(A),Q);
    set(gca,'FontSize',16)
    xlabel('\fontsize{16}Re{\it\alpha}')
    ylabel('\fontsize{16}Im{\it\alpha}')
    view(30,50);
    text(0,0,1,'\fontsize{16}(d)');
%     title('\fontsize{12}Q-function');
    subplot(245);mesh(real(A),imag(A),real(W))
    set(gca,'FontSize',16)
    xlabel('\fontsize{16}Re{\it\alpha}')
    ylabel('\fontsize{16}Im{\it\alpha}')
    view(30,50);
    text(0,0,1,'\fontsize{16}(e)');
%     title('\fontsize{12}Wigner function');
end
